HEAD ======= >>>>>>> e39bf25c59cba56a453798669eb224ec741eec47
library(tidyverse)
#library(fishmethods)
Test the Fmsy equation: Fmsy = r * (1 / (1+p)) - from the catchmsy function in R
fmsy = 0.2 * (1/(1+0.2))
#{.tabset}
##Setup
library(tidyverse)
##Fmsy Equation
Test the Fmsy equation: Fmsy = r * (1 / (1+p)) - from the catchmsy function in R
fmsy = 0.2 * (1/(1+0.2))
##Goals and Model
>>>>>>> e39bf25c59cba56a453798669eb224ec741eec47Initial Stock Assessment Test Goal: Simulate error around estimates of intial biomass
Write the function to include error in the biomass estimates. In this experiment the b_err term represents the estimated biomass with error. The biomass is randomly selected from a distribution where the mean is the actual biomass from the model and the standard deviation is 5% of the actual biomass from the model.
The basic function that we’ll use in the rest of these tests:
#Write the function
sim_stockerror <- function(b, r, r_s, error, p, k, years){
results <- data.frame(
b = rep(NA, years), c = rep(NA, years),
year = 1:years, r = rep(NA, years), b_err = rep(NA, years)) #Setup the results dataframe
#Set the initial result for the outputs in year 1
results$b[1] = b
results$b_err[1] = rnorm(1, mean = results$b[1], sd = (error*results$b[1])) #simulate error in assessment
results$c[1] = (r*(1/(1+p))) * results$b_err[1] #using f=Fmsy and b=estimate from initial stock assessment
results$r[1] = r
#Loop the model over the specified number of years
for (t in 2:years) {
results$b_err[t] = rnorm(1, mean = results$b[t-1], sd = (error*results$b[t-1]))
results$c[t] = (r*(1/(1+p))) * results$b_err[t]
results$b[t] = results$b[t-1] + (results$r[t-1] / p)*results$b[t-1]*(1 - ((results$b[t-1]/k) ^ p))-results$c[t]
results$r[t] = results$r[1] * (1 + (r_s*(t-1)))
}
return(results)
}
#Test the results
results = sim_stockerror(b = 1000, r = 0.2, error = 0.05, p = 0.2,
years = 100, r_s = 0, k = 10000)
plot(results$year, results$b)
<<<<<<< HEAD
plot(results$year, results$b_err)
Apply the new function to a set of experiments that vary actual initial biomass. In these experiments catch is responding to underlying changes in actual biomass but with some error. - Each test will change only the standard deviation associated with the error term
Standard deviation = 0.05
=======plot(results$year, results$b_err)
##Goals Tests 1-4 Apply the new function to a set of experiments that vary actual initial biomass. In these experiments catch is responding to underlying changes in actual biomass but with some error. - Each test will change only the standard deviation associated with the error term
##Test 1 Standard deviation = 0.05
>>>>>>> e39bf25c59cba56a453798669eb224ec741eec47n_experiment = 10
#Input variables
experiment <- list(
b = seq(1000, 10000, 1000),
r = rep(0.2, n_experiment),
r_s = rep(0, n_experiment),
error = rep(0.05, n_experiment),
p = rep(0.2, n_experiment),
k = rep(10000, n_experiment),
years = rep(100, n_experiment))
#Now apply to the function
results_test1 <- pmap(experiment, sim_stockerror)
#Plot results to see variation:
for(i in 1:n_experiment){
plot(results_test1[[i]]$b)
plot(results_test1[[i]]$b_err)
}
<<<<<<< HEAD
Results: Looks like regardless of initial biomass, with a error of only 5% biomass stays around Bmsy for all 100 years
Standard deviation to 0.10
======= Results: Looks like regardless of initial biomass, with a error of only 5% biomass stays around Bmsy for all 100 years
##Test 2 Standard deviation to 0.10
>>>>>>> e39bf25c59cba56a453798669eb224ec741eec47n_experiment = 10
#Input variables
experiment <- list(
b = seq(1000, 10000, 1000),
r = rep(0.2, n_experiment),
r_s = rep(0, n_experiment),
error = rep(0.1, n_experiment),
p = rep(0.2, n_experiment),
k = rep(10000, n_experiment),
years = rep(100, n_experiment))
#Now apply to the function
results_test2 <- pmap(experiment, sim_stockerror)
#Plot results to see variation:
for(i in 1:n_experiment){
plot(results_test2[[i]]$b)
plot(results_test2[[i]]$b_err)
}
<<<<<<< HEAD
Results: About the same, regardless of initial biomass, with a 10% error biomass stays around Bmsy for all 100 years
Standard deviation to 0.15
======= Results: About the same, regardless of initial biomass, with a 10% error biomass stays around Bmsy for all 100 years
##Test 3 Standard deviation to 0.15
>>>>>>> e39bf25c59cba56a453798669eb224ec741eec47n_experiment = 10
#Input variables
experiment <- list(
b = seq(1000, 10000, 1000),
r = rep(0.2, n_experiment),
r_s = rep(0, n_experiment),
error = rep(0.15, n_experiment),
p = rep(0.2, n_experiment),
k = rep(10000, n_experiment),
years = rep(100, n_experiment))
#Now apply to the function
results_test3 <- pmap(experiment, sim_stockerror)
#Plot results to see variation:
for(i in 1:n_experiment){
plot(results_test3[[i]]$b)
plot(results_test3[[i]]$b_err)
}
<<<<<<< HEAD
Results: About the same, regardless of initial biomass, with a 15% error biomass stays around Bmsy for all 100 years
Standard deviation to 0.2
======= Results: About the same, regardless of initial biomass, with a 15% error biomass stays around Bmsy for all 100 years
##Test 4 Standard deviation to 0.2
>>>>>>> e39bf25c59cba56a453798669eb224ec741eec47n_experiment = 10
#Input variables
experiment <- list(
b = seq(1000, 10000, 1000),
r = rep(0.2, n_experiment),
r_s = rep(0, n_experiment),
error = rep(0.2, n_experiment),
p = rep(0.2, n_experiment),
k = rep(10000, n_experiment),
years = rep(100, n_experiment))
#Now apply to the function
results_test4 <- pmap(experiment, sim_stockerror)
#Plot results to see variation:
for(i in 1:n_experiment){
plot(results_test4[[i]]$b)
plot(results_test4[[i]]$b_err)
}
<<<<<<< HEAD
Results: About the same, regardless of initial biomass, with a 20% error biomass stays around Bmsy for all 100 years
Use a wider range of error terms and systematically add in other variations: biomass, r, and r_s
Varying biomass with a wider range of error terms
======= Results: About the same, regardless of initial biomass, with a 20% error biomass stays around Bmsy for all 100 years
##Goals Tests 5-7 Use a wider range of error terms and systematically add in other variations: biomass, r, and r_s
##Test 5 Varying biomass with a wider range of error terms
>>>>>>> e39bf25c59cba56a453798669eb224ec741eec47First create the correct list of error terms to input in the model:
#First try creating a dataframe with the sequences we want to repeat:
input <- seq(0.1, 0.55, 0.05)
#Create a list for the output
error_input = data.frame(error = rep(NA, 100))
#Loop the repetition of each input number
for(i in 1:100){
if(i <= 10){
error_input$error[i] = input[1]
}
if(i > 10){
error_input$error[i] = input[2]
}
if(i > 20){
error_input$error[i] = input[3]
}
if(i>30){
error_input$error[i] = input[4]
}
if(i>40){
error_input$error[i] = input[5]
}
if(i>50){
error_input$error[i] = input[6]
}
if(i>60){
error_input$error[i] = input[7]
}
if(i>70){
error_input$error[i] = input[8]
}
if(i>80){
error_input$error[i] = input[9]
}
if(i>90){
error_input$error[i] = input[10]
}
}
#This works but it seems very inefficient - is there a better/faster way to get this result?
Run the experiment:
n_experiment = 100
#Input variables
experiment <- list(
b = rep(seq(1000, 10000, 1000), 10),
r = rep(0.2, n_experiment),
r_s = rep(0, n_experiment),
error = error_input$error,
p = rep(0.2, n_experiment),
k = rep(10000, n_experiment),
years = rep(100, n_experiment))
#Now apply to the function
results_test5 <- pmap(experiment, sim_stockerror)
<<<<<<< HEAD